library("tidyverse")
library("ggpubr")
library("zoo")
setwd("/mnt/LocalData/behaviour/aDN/aDN_behaviour")
gg_color_hue <- function(n) {
hues = seq(15, 375, length = n + 1)
hcl(h = hues, l = 65, c = 100)[1:n]
}
gg_color_hue(2)
[1] "#F8766D" "#00BFC4"
genotypes <- read_tsv("../2019_03_06_Courtship/genotype.tsv",col_names = TRUE)
genotypes
# indices_list <- list.files("raw data/") %>% str_subset("_Indices.csv")
indices_list <- list.files("../2019_03_06_Courtship/",recursive = TRUE) %>% str_subset("_Indices.csv") %>% str_subset("Male")
all_indices <- tibble()
for (indices_file in indices_list) {
video_name <- indices_file %>% str_remove("/.*")
temp <- read_csv(paste0("../2019_03_06_Courtship/",indices_file))
temp <- temp %>%
mutate(video = video_name)
all_indices <- bind_rows(all_indices,temp)
}
all_indices
all_male_indices <- left_join(x = genotypes,y = all_indices,by=c("video"="video","fly_id"="FlyId"))
ggplot(all_male_indices,aes(x=genotype,y=CourtshipIndexWithFacing)) +
geom_boxplot()
p1 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipIndex)) +
geom_boxplot()
p2 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipIndexWithFacing)) +
geom_boxplot()
p3 <- ggplot(all_male_indices,aes(x=genotype,y=TotalCCI)) +
geom_boxplot()
p4 <- ggplot(all_male_indices,aes(x=genotype,y=TotalCCIwFacing)) +
geom_boxplot()
ggarrange(plotlist = list(p1,p2,p3,p4),
# labels = c("ApproachingIndex",
# "ContactIndex",
# "EncirclingIndex",
# "FacingIndex",
# "TurningIndex",
# "WingIndex"),
# hjust = 1,
ncol = 4,
nrow = 1)
p1 <- ggplot(all_male_indices,aes(x=genotype,y=ApproachingIndex)) +
geom_boxplot()
p2 <- ggplot(all_male_indices,aes(x=genotype,y=ContactIndex)) +
geom_boxplot()
p3 <- ggplot(all_male_indices,aes(x=genotype,y=EncirclingIndex)) +
geom_boxplot()
p4 <- ggplot(all_male_indices,aes(x=genotype,y=FacingIndex)) +
geom_boxplot()
p5 <- ggplot(all_male_indices,aes(x=genotype,y=TurningIndex)) +
geom_boxplot()
p6 <- ggplot(all_male_indices,aes(x=genotype,y=WingIndex)) +
geom_boxplot()
ggarrange(plotlist = list(p1,p2,p3,p4,p5,p6),
# labels = c("ApproachingIndex",
# "ContactIndex",
# "EncirclingIndex",
# "FacingIndex",
# "TurningIndex",
# "WingIndex"),
# hjust = 1,
ncol = 3,
nrow = 2)
p1 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipInitiation)) +
geom_boxplot()
p2 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipTermination)) +
geom_boxplot()
p3 <- ggplot(all_male_indices,aes(x=genotype,y=CourtshipDuration)) +
geom_boxplot()
ggarrange(plotlist = list(p1,p2,p3),
ncol = 3,
nrow = 1)
df <- all_male_indices %>%
select(genotype,CourtshipTermination) %>%
group_by(genotype) %>%
mutate(len=length(CourtshipTermination))
#ggplot(df,aes(x=CourtshipTermination,color=genotype)) + geom_step(aes(len=len,y=..y.. * len),stat="ecdf")
ggplot(df,aes(x=CourtshipTermination,color=genotype)) + geom_step(aes(y=..y..),stat="ecdf")
df <- all_male_indices %>%
select(genotype,CourtshipDuration) %>%
group_by(genotype) %>%
mutate(len=length(CourtshipDuration))
#ggplot(df,aes(x=CourtshipDuration,color=genotype)) + geom_step(aes(len=len,y=..y.. * len),stat="ecdf")
ggplot(df,aes(x=CourtshipDuration,color=genotype)) + geom_step(aes(y=..y..),stat="ecdf")
p1 <- ggplot(all_male_indices,aes(x=genotype,y=ApproachingBoutLength)) +
geom_boxplot()
p2 <- ggplot(all_male_indices,aes(x=genotype,y=ContactBoutLength)) +
geom_boxplot()
p3 <- ggplot(all_male_indices,aes(x=genotype,y=EncirclingBoutLength)) +
geom_boxplot()
p4 <- ggplot(all_male_indices,aes(x=genotype,y=FacingBoutLength)) +
geom_boxplot()
p5 <- ggplot(all_male_indices,aes(x=genotype,y=TurningBoutLength)) +
geom_boxplot()
p6 <- ggplot(all_male_indices,aes(x=genotype,y=WingBoutLength)) +
geom_boxplot()
p7 <- ggplot(all_male_indices,aes(x=genotype,y=ApproachingBoutInterval)) +
geom_boxplot()
p8 <- ggplot(all_male_indices,aes(x=genotype,y=ContactBoutInterval)) +
geom_boxplot()
p9 <- ggplot(all_male_indices,aes(x=genotype,y=EncirclingBoutInterval)) +
geom_boxplot()
p10 <- ggplot(all_male_indices,aes(x=genotype,y=FacingBoutInterval)) +
geom_boxplot()
p11 <- ggplot(all_male_indices,aes(x=genotype,y=TurningBoutInterval)) +
geom_boxplot()
p12 <- ggplot(all_male_indices,aes(x=genotype,y=WingBoutInterval)) +
geom_boxplot()
ggarrange(plotlist = list(p1,p2,p3,p4,p5,p6,
p7,p8,p9,p10,p11,p12),
# labels = c("ApproachingIndex",
# "ContactIndex",
# "EncirclingIndex",
# "FacingIndex",
# "TurningIndex",
# "WingIndex"),
# hjust = 1,
ncol = 6,
nrow = 2)
#all_rawdata <- full_join(x = all_rawdata, y = genotypes, by = c("FileName"="video","Id"="fly_id"))
all_rawdata <- all_rawdata %>%
full_join(x = all_rawdata, y = genotypes, by = c("FileName"="video","Id"="fly_id")) %>%
replace_na(list(genotype = "CS female"))
# sum(is.na(all_rawdata$genotype))
# summary(all_rawdata$genotype)
# unique(all_rawdata$genotype)
all_rawdata
# # Extract density data to average
# p <- ggplot_build(test_plot)
# ggplot(as.data.frame(p$data[[1]]), aes(x,y)) + geom_line()
# ggplot(as.data.frame(p$data[[2]]), aes(x,y)) + geom_line()
all_rawdata <- all_rawdata %>%
unite("unique_fly",FileName,Id, remove = FALSE) %>%
group_by(unique_fly) %>%
mutate(
Multitasking = (Approaching + Encircling + Contact + Turning + WingGesture),
MultitaskingWithFacing = (Approaching + Encircling + Facing + Contact + Turning + WingGesture),
Courtship = ifelse(Multitasking>=1, 1, 0),
CourtshipWithFacing = ifelse(MultitaskingWithFacing>=1, 1, 0),
MultitaskingWithCopulation = (Approaching + Encircling + Contact + Turning + WingGesture + Copulation),
MultitaskingWithCopulationWithFacing = (Approaching + Encircling + Facing + Contact + Turning + WingGesture + Copulation),
CourtshipAndCopulation = ifelse(MultitaskingWithCopulation>=1, 1, 0),
CourtshipAndCopulationWthFacing = ifelse(MultitaskingWithCopulationWithFacing>=1, 1, 0),
SmoothedCourtship = ifelse((rollmean(Courtship, 150, fill = c(0,0,0), align = c("left")))>0.5, 1, 0),
SmoothedCopulation = ifelse((rollmean(Copulation, 1250, fill = c(0,0,0), align = c("center")))>0.5, 1, 0),
SmoothedDistToOther = ifelse((rollmean(ifelse(dist_to_other__mm > 2, 1, 0), 250, fill = c(1,1,NA), align = c("center")))>0.5, 1, 0)
)
test_stat_tibble0 <- all_rawdata %>%
filter(genotype != "CS female") %>%
group_by(genotype) %>%
filter(dist_to_other__mm > 2) %>%
#filter(max_wing_ang__rad > (35*pi/180)) %>%
unite("individual", FileName:Id, remove = FALSE) %>%
group_by(individual) %>%
summarise(median_facing_angle = median(facing_angle__rad),
genotype = unique(genotype))
aov0 <- aov(median_facing_angle~genotype,data = test_stat_tibble0)
summary(aov0)
Df Sum Sq Mean Sq F value Pr(>F)
genotype 3 1.036 0.3453 2.736 0.0468 *
Residuals 114 14.388 0.1262
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
TukeyHSD(aov0)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = median_facing_angle ~ genotype, data = test_stat_tibble0)
$genotype
diff lwr upr p adj
B-A -0.07051500 -0.3361762470 0.1951462 0.8999385
C-A 0.04970748 -0.1803176760 0.2797326 0.9426712
D-A 0.19276677 -0.0389146073 0.4244482 0.1380793
C-B 0.12022248 -0.1408020849 0.3812470 0.6275505
D-B 0.26328177 0.0007965104 0.5257670 0.0490063
D-C 0.14305929 -0.0832904123 0.3694090 0.3562046
test_stat_tibble0 <- all_rawdata %>%
filter(genotype != "CS female") %>%
filter(dist_to_other__mm > 2) %>%
filter(max_wing_ang__rad > (35*pi/180)) %>%
group_by(unique_fly) %>%
summarise(counts = 100*sum(facing_angle__rad > pi/8)/length(Frame),
genotype = unique(genotype))
aov0 <- aov(counts~genotype,data = test_stat_tibble0)
#summary(aov0)
TukeyHSD(aov0)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = counts ~ genotype, data = test_stat_tibble0)
$genotype
diff lwr upr p adj
B-A 6.6367127 -2.5650354 15.838461 0.2421768
C-A 7.3146600 -0.6527555 15.282076 0.0840396
D-A 16.6152204 8.5904380 24.640003 0.0000022
C-B 0.6779474 -8.3631993 9.719094 0.9973376
D-B 9.9785077 0.8867666 19.070249 0.0254776
D-C 9.3005603 1.4604520 17.140669 0.0131217
pairwise.t.test(test_stat_tibble0$counts, test_stat_tibble0$genotype)
Pairwise comparisons using t tests with pooled SD
data: test_stat_tibble0$counts and test_stat_tibble0$genotype
A B C
B 0.125 - -
C 0.055 0.845 -
D 2.2e-06 0.020 0.012
P value adjustment method: holm
test_stat_tibble0 <- all_rawdata %>%
filter(genotype != "CS female") %>%
filter(dist_to_other__mm > 2) %>%
filter(dist_to_other__mm < 10) %>%
filter(max_wing_ang__rad > (35*pi/180)) %>%
group_by(unique_fly) %>%
summarise(counts = 100*sum(facing_angle__rad > pi/8)/length(Frame),
genotype = unique(genotype))
aov0 <- aov(counts~genotype,data = test_stat_tibble0)
#summary(aov0)
TukeyHSD(aov0)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = counts ~ genotype, data = test_stat_tibble0)
$genotype
diff lwr upr p adj
B-A 4.2169183 -3.2364908 11.670328 0.4558215
C-A 4.8104624 -1.6431383 11.264063 0.2159140
D-A 11.2681432 4.7680753 17.768211 0.0000889
C-B 0.5935441 -6.7297782 7.916866 0.9966463
D-B 7.0512248 -0.3130788 14.415528 0.0658244
D-C 6.4576807 0.1071987 12.808163 0.0446836
pairwise.t.test(test_stat_tibble0$counts, test_stat_tibble0$genotype)
Pairwise comparisons using t tests with pooled SD
data: test_stat_tibble0$counts and test_stat_tibble0$genotype
A B C
B 0.286 - -
C 0.163 0.833 -
D 9.1e-05 0.056 0.046
P value adjustment method: holm
test_stat_tibble0 <- all_rawdata %>%
filter(genotype != "CS female") %>%
filter(facing_angle__rad <= pi/8) %>%
filter(max_wing_ang__rad > (25*pi/180)) %>%
group_by(unique_fly) %>%
summarise(counts = 100*sum(dist_to_other__mm > 2 & dist_to_other__mm <= 3)/length(Frame),
genotype = unique(genotype))
aov0 <- aov(counts~genotype,data = test_stat_tibble0)
#summary(aov0)
TukeyHSD(aov0)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = counts ~ genotype, data = test_stat_tibble0)
$genotype
diff lwr upr p adj
B-A 2.7255382 -8.488600 13.93968 0.9209524
C-A 2.5367905 -7.173071 12.24665 0.9040570
D-A 21.8215793 12.041805 31.60135 0.0000003
C-B -0.1887477 -11.207161 10.82967 0.9999678
D-B 19.0960411 8.015968 30.17611 0.0000986
D-C 19.2847888 9.730076 28.83950 0.0000040
test_stat_tibble <- all_rawdata %>%
filter(genotype != "CS female") %>%
#filter(genotype != "A") %>%
group_by(genotype) %>%
group_by(unique_fly) %>%
summarise(genotype = unique(genotype),
both_wing_index = 100*sum(min_wing_ang__rad[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
)
]>(35*pi/180),
na.rm = TRUE)/
length(Frame[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
)
]
))
aov1 <- aov(both_wing_index~genotype, data = test_stat_tibble)
TukeyHSD(aov1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = both_wing_index ~ genotype, data = test_stat_tibble)
$genotype
diff lwr upr p adj
B-A -0.4779248 -2.0714700 1.115620 0.8624926
C-A -0.1148557 -1.4946410 1.264930 0.9963700
D-A 0.8787432 -0.5109769 2.268463 0.3557950
C-B 0.3630691 -1.2026635 1.928802 0.9304240
D-B 1.3566679 -0.2178265 2.931162 0.1170595
D-C 0.9935989 -0.3641396 2.351337 0.2304436
test_stat_tibble <- all_rawdata %>%
filter(genotype != "CS female") %>%
#filter(genotype != "A") %>%
group_by(genotype) %>%
group_by(unique_fly) %>%
summarise(genotype = unique(genotype),
both_wing_index = 100*sum(min_wing_ang__rad[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
)
]>(15*pi/180),
na.rm = TRUE)/
length(Frame[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
)
]
))
aov1 <- aov(both_wing_index~genotype, data = test_stat_tibble)
TukeyHSD(aov1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = both_wing_index ~ genotype, data = test_stat_tibble)
$genotype
diff lwr upr p adj
B-A 1.3374241 -4.537783 7.212632 0.9338504
C-A 1.9198669 -3.167234 7.006968 0.7588352
D-A 10.3410421 5.217313 15.464771 0.0000040
C-B 0.5824427 -5.190223 6.355108 0.9935904
D-B 9.0036180 3.198649 14.808587 0.0005495
D-C 8.4211753 3.415359 13.426992 0.0001503
all_rawdata %>%
filter(genotype != "CS female") %>%
group_by(genotype) %>%
group_by(unique_fly) %>%
summarise(genotype = unique(genotype),
mean_min_wing = mean(min_wing_ang__rad[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
)
])
) %>%
ggplot(aes(x=genotype,y=mean_min_wing)) +
geom_boxplot() +
ylim(0,.5)
all_rawdata %>%
filter(genotype != "CS female") %>%
group_by(genotype) %>%
group_by(unique_fly) %>%
summarise(genotype = unique(genotype),
min_wing_gt_15 = sum(mean(min_wing_ang__rad[which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
)
])>15*pi/180)
) %>%
summarise(prop_min_wing_gt_15 = 100*sum(min_wing_gt_15>0,na.rm = TRUE)/length(min_wing_gt_15))
all_rawdata %>%
filter(dist_to_other__mm > 2) %>%
filter(dist_to_other__mm < 10) %>%
filter(max_wing_ang__rad > (25*pi/180)) %>%
#filter(genotype == "D") %>%
ggplot(aes(x=genotype, y=wing_l_len__px))+
geom_boxplot()
Warning messages:
1: Unknown or uninitialised column: 'bin_max_wing'.
2: Unknown or uninitialised column: 'bin_max_wing'.
3: Unknown or uninitialised column: 'bin_max_wing'.
4: Unknown or uninitialised column: 'bin_max_wing'.
5: Unknown or uninitialised column: 'bin_max_wing'.
6: Unknown or uninitialised column: 'bin_max_wing'.
7: Unknown or uninitialised column: 'bin_max_wing'.
8: Unknown or uninitialised column: 'bin_max_wing'.
temp <- all_rawdata %>%
filter(genotype!="CS female") %>%
filter(dist_to_other__mm > 2)
uniq_fly <- unique(temp$unique_fly)
facing_start_and_end_wing <- tibble()
for (fly in uniq_fly) {
temp2 <- filter(temp, unique_fly == fly)
bouts <- rle(temp2$WingGesture)
starts <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
ends <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
ends[i] = sum(bouts$lengths[1:(i*2)])
}
genotype = rep(unique(temp2$genotype),length(which(bouts$values==1)))
uni_id = rep(unique(temp2$unique_fly),length(which(bouts$values==1)))
temp_tibble <- tibble(unique_fly = uni_id,
genotype = genotype,
facing_at_start = temp2$facing_angle__rad[starts],
facing_at_end = temp2$facing_angle__rad[ends]
)
facing_start_and_end_wing <- bind_rows(facing_start_and_end_wing,temp_tibble)
}
facing_start_and_end_wing
temp <- all_rawdata %>%
filter(genotype!="CS female") %>%
filter(dist_to_other__mm > 2)
uniq_fly <- unique(temp$unique_fly)
facing_start_and_end_wing <- tibble()
for (fly in uniq_fly) {
temp2 <- filter(temp, unique_fly == fly)
bouts <- rle(temp2$WingGesture)
starts <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
ends <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
ends[i] = sum(bouts$lengths[1:(i*2)])
}
genotype = unique(temp2$genotype)
uni_id = unique(temp2$unique_fly)
temp_tibble <- tibble(unique_fly = uni_id,
genotype = genotype,
facing_at_start = mean(temp2$facing_angle__rad[starts]),
facing_at_end = mean(temp2$facing_angle__rad[ends])
)
facing_start_and_end_wing <- bind_rows(facing_start_and_end_wing,temp_tibble)
}
facing_start_and_end_wing
temp <- all_rawdata %>%
filter(genotype!="CS female") %>%
filter(dist_to_other__mm > 2)
uniq_fly <- unique(temp$unique_fly)
facing_start_and_end_wing <- tibble()
for (fly in uniq_fly) {
temp2 <- filter(temp, unique_fly == fly)
bouts <- rle(temp2$WingGesture)
starts <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
ends <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
ends[i] = sum(bouts$lengths[1:(i*2)])
}
genotype = unique(temp2$genotype)
uni_id = unique(temp2$unique_fly)
temp_tibble <- tibble(unique_fly = uni_id,
genotype = genotype,
facing_at_start = median(temp2$facing_angle__rad[starts]),
facing_at_end = median(temp2$facing_angle__rad[ends])
)
facing_start_and_end_wing <- bind_rows(facing_start_and_end_wing,temp_tibble)
}
facing_start_and_end_wing
temp <- all_rawdata %>%
filter(genotype!="CS female") %>%
filter(dist_to_other__mm > 2)
Warning messages:
1: Unknown or uninitialised column: 'bin_max_wing'.
2: Unknown or uninitialised column: 'bin_max_wing'.
3: Unknown or uninitialised column: 'bin_max_wing'.
4: Unknown or uninitialised column: 'bin_max_wing'.
uniq_fly <- unique(temp$unique_fly)
facing_start_and_end_wing <- tibble()
for (fly in uniq_fly) {
temp2 <- filter(temp, unique_fly == fly)
#bouts <- rle(temp2$WingGesture)
temp2$bin_max_wing <- ifelse(temp2$max_wing_ang__rad >= (35*pi/180),1,0)
temp2$bin_max_wing <- temp2$bin_max_wing %>% replace_na(0)
bouts <- rle(temp2$bin_max_wing)
starts <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
ends <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
ends[i] = sum(bouts$lengths[1:(i*2)])
}
genotype = unique(temp2$genotype)
uni_id = unique(temp2$unique_fly)
temp_tibble <- tibble(unique_fly = uni_id,
genotype = genotype,
facing_at_start = median(temp2$facing_angle__rad[starts]),
facing_at_end = median(temp2$facing_angle__rad[ends])
)
facing_start_and_end_wing <- bind_rows(facing_start_and_end_wing,temp_tibble)
}
facing_start_and_end_wing
test_stat_tibble0 <- all_rawdata %>%
filter(genotype!="CS female") %>%
group_by(unique_fly) %>%
filter(dist_to_other__mm > 2) %>%
filter(max_wing_ang__rad > (35*pi/180)) %>%
summarise(genotype = unique(genotype),
facing = median(facing_angle__rad))
aov0 <- aov(facing~genotype,data = test_stat_tibble0)
#summary(aov0)
TukeyHSD(aov0)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = facing ~ genotype, data = test_stat_tibble0)
$genotype
diff lwr upr p adj
B-A -0.064278836 -0.2941883 0.1656306 0.8852971
C-A -0.017473865 -0.2165430 0.1815953 0.9957523
D-A -0.012774271 -0.2132767 0.1877282 0.9983603
C-B 0.046804972 -0.1790918 0.2727017 0.9489797
D-B 0.051504566 -0.1756563 0.2786654 0.9345696
D-C 0.004699594 -0.1911887 0.2005879 0.9999116
test <- all_rawdata %>%
filter(FileName == "Megan-2019_03_06_Courtship-DsxVglutTNT_Male_1234_2") %>%
filter(Id == 36)
temp2 <- test
Warning messages:
1: Unknown or uninitialised column: 'bin_max_wing'.
2: Unknown or uninitialised column: 'bin_max_wing'.
3: Unknown or uninitialised column: 'bin_max_wing'.
4: Unknown or uninitialised column: 'bin_max_wing'.
#bouts <- rle(temp2$WingGesture)
temp2$bin_max_wing <- ifelse(temp2$max_wing_ang__rad > (35*pi/180),1,0)
temp2$bin_max_wing <- temp2$bin_max_wing %>% replace_na(0)
bouts <- rle(temp2$bin_max_wing)
starts <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
ends <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
ends[i] = sum(bouts$lengths[1:(i*2)])
}
genotype = rep(unique(temp2$genotype),length(which(bouts$values==1)))
uni_id = rep(unique(temp2$unique_fly),length(which(bouts$values==1)))
temp_tibble <- tibble(unique_fly = uni_id,
genotype = genotype,
facing_at_start = temp2$facing_angle__rad[starts],
facing_at_end = temp2$facing_angle__rad[ends]
)
temp_tibble
starts = for (i in seq(1,length(which(bouts$values==1)),1)) {
starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
starts
NULL
test$bin_max_wing <- ifelse(test$max_wing_ang__rad > (35*pi/180),1,0)
test$bin_max_wing <- test$bin_max_wing %>% replace_na(0)
#test$bin_max_wing
bouts <- rle(test$bin_max_wing)
bouts
Run Length Encoding
lengths: int [1:393] 675 3 2 14 2 22 212 15 15 2 ...
values : num [1:393] 0 1 0 1 0 1 0 1 0 1 ...
bouts <- rle(test$WingGesture)
bouts
Run Length Encoding
lengths: int [1:1361] 666 5 3 45 198 1 2 5 4 17 ...
values : num [1:1361] 0 1 0 1 0 1 0 1 0 1 ...
sum(test$bin_max_wing[which.max(test$SmoothedCourtship):which.max(test$SmoothedCopulation)]==test$WingGesture[which.max(test$SmoothedCourtship):which.max(test$SmoothedCopulation)])
[1] 1230
sum(test$bin_max_wing[which.max(test$SmoothedCourtship):which.max(test$SmoothedCopulation)]!=test$WingGesture[which.max(test$SmoothedCourtship):which.max(test$SmoothedCopulation)])
[1] 35
gg_color_hue(4)
[1] "#F8766D" "#7CAE00" "#00BFC4" "#C77CFF"
max(test$Frame)
[1] 5974
temp2 <- all_rawdata %>%
filter(FileName == "Megan-2019_03_06_Courtship-DsxVglutTNT_Male_123_3") %>%
filter(Id == 1) %>%
filter(dist_to_other__mm > 2)
#bouts <- rle(temp2$WingGesture)
temp2$bin_max_wing <- ifelse(temp2$max_wing_ang__rad >= (35*pi/180),1,0)
temp2$bin_max_wing <- temp2$bin_max_wing %>% replace_na(0)
bouts <- rle(test$bin_max_wing)
starts <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
starts[i] = sum(bouts$lengths[1:(i*2)-1])+1
}
ends <- as.numeric()
for (i in seq(1,length(which(bouts$values==1)),1)) {
ends[i] = sum(bouts$lengths[1:(i*2)])
}
genotype = unique(temp2$genotype)
uni_id = unique(temp2$unique_fly)
temp_tibble <- tibble(unique_fly = uni_id,
genotype = genotype,
facing_at_start = median(temp2$facing_angle__rad[starts]),
facing_at_end = median(temp2$facing_angle__rad[ends])
)
temp_tibble
bouts <- rle(test$bin_max_wing)
bouts
Run Length Encoding
lengths: int [1:198] 631 3 57 5 363 34 357 1 23 1 ...
values : num [1:198] 0 1 0 1 0 1 0 1 0 1 ...
If i have to for loop over the groups of a tibble, the below line might be usefull. dim(unique(all_rawdata[all_rawdata %>% group_by(genotype) %>% group_vars()]))[1]
courtship_window <- function(input,wind=600,...){
temp_tibble <- summarise(input, start_of_courtship = which.max(SmoothedCourtship),
end_of_courtship = ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
)
)
return(temp_tibble$start_of_courtship:temp_tibble$end_of_courtship)
#return(c(temp_tibble$start_of_courtship,temp_tibble$end_of_courtship))
#return(temp_tibble)
}
courtship_window <- function(input,wind=600,...){
temp_tibble <- summarise(input, start_of_courtship = which.max(SmoothedCourtship),
end_of_courtship = ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
)
)
#return(temp_tibble$start_of_courtship:temp_tibble$end_of_courtship)
return(c(temp_tibble$start_of_courtship,temp_tibble$end_of_courtship))
#return(temp_tibble)
}
courtship_window(test, wind = 600)
temp <- all_rawdata %>%
filter(genotype!="CS female")
uniq_fly <- unique(temp$unique_fly)
courting_frames <- tibble()
for (fly in uniq_fly) {
temp2 <- temp %>% filter(unique_fly == fly)
temp2 <- slice(temp2, which.max(SmoothedCourtship):ifelse(which.max(SmoothedCopulation) > which.max(SmoothedCourtship),
ifelse(which.max(SmoothedCopulation) <= (which.max(SmoothedCourtship)+(25*wind)),
which.max(SmoothedCopulation),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
),
min((which.max(SmoothedCourtship)+(25*wind)),max(Frame))
))
courting_frames <- bind_rows(courting_frames,temp2)
}
courting_frames
calculate_indices <- function(input = .,feature,jaaba=TRUE,thresh=NULL,wind=600){
if (feature %in% names(input)) {
var_ind = paste0(feature,"_index")
if (jaaba) {
temp_tibble <- summarise(input, genotype = unique(genotype),
!!var_ind := 100*sum(input[[feature]][which.max(input[["SmoothedCourtship"]]):ifelse(which.max(input[["SmoothedCopulation"]]) > which.max(input[["SmoothedCourtship"]]),
ifelse(which.max(input[["SmoothedCopulation"]]) <= (which.max(input[["SmoothedCourtship"]])+(25*wind)),
which.max(input[["SmoothedCopulation"]]),
min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
),
min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
)])/
length(input[["Frame"]][which.max(input[["SmoothedCourtship"]]):ifelse(which.max(input[["SmoothedCopulation"]]) > which.max(input[["SmoothedCourtship"]]),
ifelse(which.max(input[["SmoothedCopulation"]]) <= (which.max(input[["SmoothedCourtship"]])+(25*wind)),
which.max(input[["SmoothedCopulation"]]),
min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
),
min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
)])
)
} else {
temp_tibble <- summarise(input, genotype = unique(genotype),
!!var_ind := 100*sum(input[[feature]][which.max(input[["SmoothedCourtship"]]):ifelse(which.max(input[["SmoothedCopulation"]]) > which.max(input[["SmoothedCourtship"]]),
ifelse(which.max(input[["SmoothedCopulation"]]) <= (which.max(input[["SmoothedCourtship"]])+(25*wind)),
which.max(input[["SmoothedCopulation"]]),
min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
),
min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
)]>thresh,na.rm = TRUE)/
length(input[["Frame"]][which.max(input[["SmoothedCourtship"]]):ifelse(which.max(input[["SmoothedCopulation"]]) > which.max(input[["SmoothedCourtship"]]),
ifelse(which.max(input[["SmoothedCopulation"]]) <= (which.max(input[["SmoothedCourtship"]])+(25*wind)),
which.max(input[["SmoothedCopulation"]]),
min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
),
min((which.max(input[["SmoothedCourtship"]])+(25*wind)),max(input[["Frame"]]))
)])
)
}
return(temp_tibble)
} else {
message(paste0(feature, " does not exist in table"))
}
}
calculate_indices <- function(input,feature,jaaba=TRUE,thresh=NULL){
if (feature %in% names(input)) {
var_ind = paste0(feature,"_index")
if (jaaba) {
temp_tibble <- summarise(input, genotype = unique(genotype),
!!var_ind := 100*sum(input[[feature]][courtship_window(input)])/
length(input[["Frame"]][courtship_window(input)])
)
} else {
temp_tibble <- summarise(input, genotype = unique(genotype),
!!var_ind := 100*sum(input[[feature]][courtship_window(input)]>thresh,na.rm = TRUE)/
length(input[["Frame"]][courtship_window(input)])
)
}
return(temp_tibble)
} else {
message(paste0(feature, " does not exist in table"))
}
}
all_rawdata %>%
filter(genotype != "CS female") %>%
filter(genotype == "D") %>%
group_by(genotype) %>%
group_by(unique_fly) %>%
calculate_indices(feature = "max_wing_ang__rad",jaaba = FALSE,thresh = (25*pi/180))
test <- all_rawdata %>%
filter(FileName == "Megan-2019_03_06_Courtship-DsxVglutTNT_Male_1234_2") %>%
filter(Id == 36)
calculate_indices(input = test,feature = "WingGesture",jaaba = TRUE)
calculate_indices(input = test,feature = "max_wing_ang__rad",jaaba = FALSE,thresh = (25*pi/180))
calculate_indices(input = test,feature = "max_wing_angle__rad",jaaba = FALSE,thresh = (25*pi/180))